Quantum Monte Carlo study of spin-polarized deuterium 
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The ground state properties of spin-polarized deuterium (D4-) at zero temperature are obtained 
by means of the diffusion Monte Carlo calculations within the fixed-node approximation. Three D\, 
species have been investigated (D4.1, DJ,2, D^), corresponding respectively to one, two and three 
equally occupied nuclear spin states. Influence of the backflow correlations on the ground state 
energy of the systems is explored. The equilibrium densities for DJ,2 and D4.3 liquids are obtained 
and compared with ones obtained in previous approximate prediction. The density and the pressure 
at which the gas-liquid phase transition occurs at T=Q is obtained for DJ,i. 



I. INTRODUCTION 

As the simplest element in the nature, hydrogen, has 
been investigated theoretically and experimentally from 
the very first beginnings of the quantum theory appli- 
cations. Since the hydrogen appears in the three iso- 
topic forms (hydrogen, deuterium and tritium), it of- 
fers even more interesting scientific investigation possi- 
bilities. Very strong interest of the scientific commu- 
nity for the study of the electron spin-polarised hydro- 
gen (H4.) and its isotopes, spin-polarised deuterium (DJ.) 
and spin-polarised tritium (T\f), begins after Kolos and 
Wolniewicz (KW) in 1965^ calculated very precisely the 
triplet pair potential b 3 Y<+ . A reasonably huge enthusi- 
asms for the study of the electron spin-polarised hydrogen 
and its isotopes was reflected through the expectations 
to explore even more extreme quantum behaviour then 
those obtained in study of the helium isotopes 2 - - — . Such 
expectations were grounded on the extremely weak at- 
traction of the triplet pair potential b 3 T,^ through which 
the two HI (D| or T|) atoms interaction is described, 
and of course, on even lighter masses compared to the 
helium isotopes masses. 

It has been shown by Freed^ in 1980 that, within the 
Born-Oppenheimer approximation in the spin-aligned 
electronic state, hydrogen nuclei behave as effective 
bosons, as well as tritium nuclei, an even before that 
Stwalley and Nosanow 3 - in 1976 proposed HJ. as the best 
new candidate for a Bose- Einstein condensate (BEC). 
Their prediction was experimentally confirmed in 1998 
when Fried et alfi managed to form a BEC state of HJ, 
using experimental set-up that assumed a wall-free con- 
finement and a low evaporation rates. This very impor- 
tant experimental achievement was accompanied by the 
intensive theoretical predictions. Recently the ground 
state properties of HJ, have been investigated using diffu- 
sion Monte Carlo (DMC) method^. By means of the very 
precise microscopic calculations it was confirmed that HJ. 
does not form a liquid phase at zero temperature, and 
in addition the gas-solid phase transition that occurs at 
T=0 was also examined. The ground state properties 
obtained in DMC calculations are obtained using the 
triplet pair potential 6 3 £j which has been recently re- 



calculated and extended to larger interparticle distances 
by Jamieson et al. (JDW) 8 . The obtained DMC re- 
sults are in good agreement with the conclusions, pre- 
viously obtained with different variational methods by 
other authors 3 -^—, concerning the ground state HJ, gas 
phase. 

Very vivid theoretical interest was dedicated to the 
study of the ground state properties of the T| bulk. Due 
to its larger mass and the fact that T\. atoms obey Bose 
statistics, it was predicted using the variational calcula- 
tions that T4- forms a liquid ground state^i 2 -. Those 
predictions have been also recently confirmed with the 
very precise microscopic DMC calculations^ 3 -, and the 
densities at which the liquid-solid phase transition occurs 
at T=0 are also determined. As in case of spin-polarised 
hydrogen^, the JDW interatomic potential was used in 
DMC study of spin-polarised tritium. In addition to all 
previously refereed, Blume et ali^ pointed TJ, as a new 
possible BEC gas in the optical dipol trap, as a conse- 
quence of the very broad Fesbach resonance that can be 
used to tune interaction in the trap. 

Up to now, according to the best of our knowledge, 
precise microscopic DMC calculations have not been pre- 
formed yet for the spin-polarised deuterium, leaving thus 
determination of the equation of state for DJ, only at 
the variational level. Having in mind that the DJ. is a 
Fermi system with nuclear spin value 1, and the fact 
that the three DJ, species are possible depending on how 
DJ, atoms are distributed with respect to the available 
nuclear spin states (-1, 0, +1), an extensive theoretical 
study was dedicated to the ground state properties of 
the one of the most interesting system among the exist- 
ing quantum systems^—. It has been shown that when 
all three nuclear spin states are equally occupied, D4.3, 
the liquid state is the ground state of the system-^. The 
ground state of the system in which two nuclear spin 
state are equally occupied, D4.2, was finally resolved by 
Panoff and Clark^. when they used improved wave func- 
tion which included optimal Jastrow, triplet and back- 
flow correlations in their variational Monte Carlo (VMC) 
calculations. The negative energy per particle values ob- 
tained at equilibrium density po = 0.004 A" 3 and zero 
pressure revealed that D4.2 also forms a self-bound liquid 
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in the ground state. In the same work the ground state 
of the Dli system has not been fully resolved because 
of the positive varational energy per particle that was 
obtained in the case of DJ,i , even when the refinements 
to the ground state trial wave function were included in 
the description of the system. Since, they could not pre- 
dict with certainty the zero-temperature phase of D^i, 
they concluded that D4.1 may remain in the gaseous state 
down to absolute zero, but they left open possibility for 
D-li to liquefy under very slight pressure. Their results 
have been qualitatively confirmed recently by Skjetne 
and 0stgaard with the Silvera interaction potential^. In 
our recent VMC study of D|i bulk we investigated the 
influence of the interaction potential on the D4.1 equa- 
tion of state and we discussed the liquid-gas coexistence 
region 20 . Our variational calculations showed that the 
gas-liquid transition occurs at extremely low density of 
the gas (~ 10~ 5 A -3 ) and very low pressure (^0.0008 
bar). It was obvious that we need to proceed with the 
steep forward, i.e. to continue with the DMC calcula- 
tions in order to determine more precisely the gas-liquid 
transition. 



In the present work, we present our results obtained us- 
ing the diffusion Monte Carlo calculations for the three 
spin-polarised deuterium species, DJ^i, and D4.3 at 
zero-temperature. Since the Fermion systems are consid- 
ered within this study we had to deal with the well-known 
sign problem of the probe wave function, and because of 
that we employed a fixed-node (FN) approximation in 
our theoretical approach. This approximation is known 
as one of the best theoretical models for prediction of the 
ground state properties of the Fermi systems, especially 
in cases in which the nodes of the probe wave function 
are very close to the exact nodes. Interaction between 
D| is modelled using the newest JDW triplet pair po- 
tential 6 3 S„ , where a cubic spline is used to interpolate 
between JDW data. This interaction is then smoothly 
connected to the long-range behaviour stated in Ref. [2l|. 
In addition to the accurate microscopic results of the en- 
ergetic and structural properties of the D4.1, D| 2 and 
D4.3 bulk systems at T=0, we comment influence of the 
backflow correlations on the ground state energy at the 
DMC level. The obtained equilibrium densities for DJ,3 
and D4.2 liquids are compared with ones determined in 
previously approximate descriptions. The possibility of 
the gas-liquid phase transition at T=0 for D4.1 is explored 
using DMC results. 



In Sec. II, the DMC method and the trial wave func- 
tion used for importance sampling are briefly described. 
The results obtained for the three spin-polarized deu- 
terium species DJ^i, D4.2 and D4.3, equations of state as 
well as the structural properties are presented in subsec- 
tions of the Sec. III. The main conclusions of this work 
are summarised one more time in Sec. IV. 



II. METHOD 

The focus of the diffusion Monte Carlo method is solv- 
ing stochastically the Schrodinger equation written in 
imaginary time, 

-h^^- = (H-E T )nR,t) , (1) 
for an iV-particle Hamiltonian 
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The constant which in Eq. (fT]) acts as a reference en- 
ergy is E T , and R = (rj., . . . , r^) is known as a walker, 
which collectively denotes the particles positions within 
the Monte Carlo methodology. In order to continue 
with the stochastic solution, the Schrodinger equation 
has to been rewritten in terms of the mixed distribu- 
tion $(R, t) = ^(R, t)ijj(R), where a trial wave function 
tp(R) has to be introduced for the importance sampling. 
In the diffusion process the mixed distribution $(R, t) is 
represented by a set of walkers. The lowest energy eigen- 
function, not orthogonal to ip(R) survives in the limit 
t — > 00, and then the sampling of the ground state is 
effectively achieved for a TV-body bosonic system. Be- 
cause of the sign problem which exists in case of Fermi 
systems, the ^(R,t)ip(R) is not always grater or equal 
zero. To satisfy the condition ^(R, t)ip(R) > 0, it is use- 
ful to use the fixed-node approximation in which ^(R, t) 
and ip(R) have to change sign together, i.e. to share the 
same nodes. Having this in mind, the fixed node ener- 
gies obtained in the t — > 00 limit are the upper bounds to 
the exact ground state energies. Only in case of the ex- 
act nodes of the probe wave function ip(R), the obtained 
fixed node energies are the ground state energies. 

The trial wave function used in simulations is the 
Jastrow-Slater model ip(R) = ipAipJ = Tii<j f( r ij)i 
where ifjj is the part of the trial wave function which de- 
scribes the dynamical correlations induced by V(rij), and 
ipA is the antisymmetric part which introduces the statis- 
tics of Fermi particles. The Schiff-Verlet (SV)^ 2 - form 



/(r-y) = exp 



(3) 



with variational parameter b is used to model the two- 
body correlations. The antisymmetric part of the trial 
wave function ipA is modelled with a Slater determinant 
in case of DJ,i , and product of Slater determinants in case 
fo DJ,2 and DJ,3. Single-particle plane wave orbitals are 
used in the Slater determinant, ip ai (vj) — exp(i k Q< rj), 
which correspond to the exact wave function of the Fermi 
sea. 

As it is usual for the bulk system simulations with a 
finite number of particles, a size-dependence analysis has 
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also been explored. We added to our DMC results the 
standard tail corrections, i.e. corrections coming from 
the finite size of the simulation box, as well as the Fermi 
correction for the ground state energy of the system. The 
Fermi correction was obtained as the difference between 
the ground state energy per particle of a free Fermi gas 

(tlm ip^P) 3 ) an d the result obtained by summing the 
discrete contributions of wave vectors (K^-) used in our 
finite N particle simulation. As it is shown in Refl2G. 
it is enough to use 33 particles to simulate D|i bulk. 
The same values of the optimal variational parameter b 
reported in Ref[13 are used in diffusion Monte Carlo cal- 
culations of D4.1 for the investigated density range from 
0.00009 to 0.00634 A~ 3 . 

In case of DJ,2 and the density range from 0.00282 to 
0.00634 A~ 3 the optimal value of the parameter b slightly 
increases from 3.89 to 3.97 A. After we have performed 
the detailed size dependence analysis in the case of den- 
sity p = 0.00493 A" 3 , which is given in Table HI we de- 
cided to use the 66 particles in the DMC simulations. As 
it is shown in Table HI the minimum value of the Fermi 
correction is produced for the 66 particle system. 



for the 99 particles, so we decided to continue with the 
DMC calculations with the 99 atoms. 



TABLE I. Energy per particle for the D4.2 system as a function 
of the number of atoms included in the simulation (in K) 
for density p = 0.00493 A" 3 . Results are obtained with the 
VMC and the JDW interatomic potential 8 . The figures in 
parenthesis are the statistical errors. 

In case of D4,3 we have optimized 57 particle system 
for density p = 0.00352 A~ 3 and we have used that pa- 
rameter b for all densities. Namely, we have obtained the 
value b = 3.93 A which minimises energy per particle 
in case of D^ system at density p = 0.00352 A -3 , and 
Panoff and Clark have used b = 1.065a — 3.93 A for all 
considered density in case of assuming only Schiff-Verlet 
correlations 1 ^. To check is this value of the parameter 
b is really the one which optimizes the energy, we have 
performed additional optimization with the 99 particles. 
In the density range from 0.00282 to 0.00634 A" 3 the op- 
timal value of the parameter b slightly changed from 3.93 
to 4.01 A, but there was no change in the result for the 
energy per particle in the VMC calculations, compared to 
the results obtained with the parameter b used by Panoff 
and Clark. We concluded that the variational parameter 
used by Panoff and Clark is very well determined and 
we performed the DMC calculations using that param- 
eter (b = 3.93 A). Detailed size dependence analysis in 
the case of density p = 0.00282 A~ 3 is given in Table [II] 
The minimum value of the Fermi correction is obtained 



N 


E/N 


E/N 


Fermi 


E/N 




VMC 


tail 


correction 


total 


57 


0.071(10) 


-0.077 


0.041 


0.04(1) 


99 


0.101(20) 


-0.047 


0.005 


0.06(2) 


171 


0.081(10) 


-0.028 


-0.017 


0.04(1) 



TABLE II. Energy per particle for the D4.3 system as a func- 
tion of the number of atoms included in the simulation (in 
K) for density p = 0.00282 A" 3 . Results are obtained with 
the VMC and the JDW interatomic potential 8 -. The figures 
in parenthesis are the statistical errors. 



Since the results of the FN approximation depend on 
the quality of the trial wave function ip, we tried to im- 
prove the trial wave function by introducing the moment- 
dependent correlations in the trial wave function. The 
backfiow correlations we modelled in the similar way 
as Panoff and Clark dicU^, but omitting the long-range 
term, i.e. in a following way 
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E/N 


E/N 


Fermi 


E/N 




VMC 


tail 


correction 


total 


38 


0.810(10) 


-0.305 


0.079 


0.58(1) 


54 


0.841(10) 


-0.223 


0.075 


0.69(1) 


66 


0.837(10) 


-0.188 


0.01 


0.66(1) 


114 


0.830(10) 


-0.119 


-0.031 


0.68(1) 



Yj = Tj + X B eX P 



UJB 



(r,-r fe ), (4) 



where As, rs and ojb are additional variational parame- 
ters. 

At the variational Monte Carlo level the introduction 
of the backfiow correlations in the case of D4.1 could not 
improve the results because Xb — > and the lowering of 
the energy was practically zero. The variational param- 
eters that minimised the energy per particle of the D4.2 
bulk are A B =0.14, r s =2.33 A and w B =1.68 A, while 
those of the D| 3 bulk are A B =0.56, ri3=1.28 A and 
u; B =2.53 A. For all three species we have minimized the 
variational parameters of the backfloe correlations for the 
density 0.00352 A~ 3 , which is the density close to the 
minimum of the equation of state obtained with the VMC 
calculations in the case when the backfiow correlations 
are not included in the model. 

For all three spin-polarized deuterium species we used 
the DMC method accurate to the second order in the 
time-step A*2£. In order to reduce any systematic bias 
coming from the time-step or mean number of walkers 
used in simulations, we investigated carefully a possible 
dependences. 



III. RESULTS 

A. Bulk D^i equation of state 

The bulk D4.1 system is studied in the density range 
from 0.00009 to 0.00634 A~ 3 and the obtained DMC re- 
sults are plotted in Fig. [T] Since at the VMC level the 



4 



Di, 



0.45 



Dj 0.3 



0.002 0.004 

P(A" 3 ) 




40000 80000 120000 

1/p(A 3 ) 



FIG. 1. Energy per particle of D^i without backflow corre- 
lations (solid triangles) and with backflow correlations (solid 
circles) as a function of the density p. The error bars of the 
DMC energies are smaller than the size of the symbols. 

introduction of the backflow correlations in the case of 
D^-i could not improve the results, we tried with the ad- 
ditional minimizations of the backflow parameters at the 
DMC level. For the density 0.00423 A -3 , which is the 
density close to the density that defines energy minimum 
at the DMC level in the case when the backlfow cor- 
relations have not been included in the model, we have 
performed additional DMC calculations in which we have 
changed the backflow parameters Xb, Tb and ujb- In that 
way we tried to explore the quality of the nodes of our 
trial way function. This investigation showed that en- 
ergy per particle slightly lowers if the DMC calculations 
are preformed with the backflow parameters As=0.14, 
rs=2.43 A and wb=1.52 A. The lowering of the energy 
is practically negligible in the region of the small densi- 
ties and in the region of the larger densities resides up 
to 13%. Since we are dealing with the Fermi system it is 
not surprising that variational Monte Carlo was not able 
to provide the best guiding function. 

To compare our results for equilibrium density and 
energy per particle at that density with the previously 
obtained results reported in Refll8l. we interpolated our 
DMC results in the density range from 0.00282 to 0.00634 
A" 3 using the fit of the form (e = E/N) 

e{ P ) = e + B [^)\c{P-^)\ (5) 

where po and eo being respectively the equilibrium den- 
sity and energy per particle at this density. The same 
density range was investigated in Refflil and given ana- 
lytical form we used to obtain equation of state of the 
liquid spin-polarised tritium in Ref. Il3l 

In the case where the backflow correlations have not 
been included in the model the fitting resulted with the 
best set of parameters e = 0.1246(4) K, B = 1.263(7) 
K, C = 0.83(1) K, and p a = 0.004169(3) A~ 3 , and 
when the backflow correlations are included in the model 
e = 0.1086(8) K, B = 1.31(2) K, C = 0.8(1) K, 
and po = 0.00420(3) A~ 3 . It is evident that inclusion 
of the backflow correlations in the model moves equi- 



FIG. 2. Energy per particle as a function of 1/p for DJ,i- 
Dashed line represents universal expansion for Fermi gas in 
the region of very small densities. The common tangent to 
the liquid and the gas equations of state is shown as a dotted 
line. The inset shows the same results (except the common 
tangent) in logarithmic scale. 



librium density to the slightly higher density and low- 
ers eo around 13%. The comparison of our best result 
(e = 0.1086(8) K and p = 0.00420(3) A~ 3 ) with the 
best variational result of Panoff and Clark reported in 
Refill (e = 0.26(1) K and p = 0.004 A" 3 ) reveals 
also moving the equilibrium density to the slightly higher 
density and lowering of the energy per particle. As it 
is expected, lowering of the energy due to the diffusion 
process in DMC method is more significant than the low- 
ering caused by improving the probe wave function with 
backflow correlations. The conclusion would not change 
if we were refereeing to the results of the same authors 
reported in Refill (e min = 0.26 K and p mm S 0.00349 

A- 3 ) . 
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FIG. 3. Two-body radial distribution functions of D4.1. From 
bottom to top in the height of the main peak, the results 
correspond to densities 0.00423, 0.00493, 0.00563, and 00634 
A -3 . The inset shows g(r) in the case of the extremely low 
density p = 0.00052 A -3 . 

The results for energy per particle that we present 
within this work for D4.1, as well as the results obtained 
in our previous variational investigation for the JDW in- 
teraction potential^, are in the liquid-gas coexistence 
region^. In the literature liquid-gas coexistence region 
is defined as the region in which a first order gas-liquid 
phase transition is possible at absolute zero. In order 



5 



to construct a common tangent to the liquid and the 
gas equations of state, i.e. the Maxwell double-tangent, 
we first plot in Fig 12] our results as a function of 1/p, 
i.e., volume of the system. To proceed with the Maxwell 
double-tangent construction we had to include the uni- 
versal expansion for a Fermi gas in the region of very 
small densities^. Namely, it is practically impossible to 
simulate with diffusion Monte Carlo method quantum 
system in the extremely low density regime. In Figj2]the 
universal expansion for a Fermi gas in the region of very 
small densities we plotted with a dashed line, and the 
Maxwell double-tangent is shown as a dotted line. The 
presented results indicate that the gaseous state is the 
ground state of Dii , and that D^i liquefies by implying 
just a very slight pressure at T = 0. A first order gas- 
liquid transition occurs at gas density 1.48 • 10~ 5 A~ 3 
and extremely low pressure of ^0.00009 bars. With our 
former variational results we could predict the transition 
at the gas density 5.4 • 10~ 5 A -3 and slightly higher pres- 
sure of ^0.0008 bars. Both predictions of the transition 
suppose application of the extremely low pressure, that 
was already pointed in Reffl8l. 

In order to present our results more clearly we added 
the inset in Fig[2] which shows the same results in loga- 
rithmic scale. The qualitative conclusion regarding the 
liquid-gas coexistence region would not change if we have 
used a more precise description including three-body po- 
tential of interaction 2 -.. Such more sophisticated inter- 
action potential in used in the work of Blume et al. in 
the study of the ground state properties of small spin- 
polarized tritium clusters 14 . In that work Blume et al. 
showed that inclusion of the three-body potential in the 
Hamiltonian results with small increase of the ground 
state energy of the clusters. 



p=0. 00423 A' 3 




k(A"') 

FIG. 4. Static structure function of D^i- From bottom to 
top in the height of the main peak, the results correspond 
to densities 0.00423, 0.00493, 0.00563, and 00634 A" 3 . The 
inset shows S(k) in the case of the extremely low density 
p = 0.00052 A" 3 . 

The ground state structure properties of D4.1 are ob- 
tained from the two-body radial distribution function 
and from its Fourier transform, the static structure fac- 
tor. In Fig|3]and Figf4]we present the DMC results ob- 
tained at different densities using the method for the pure 
estimators^. As it is expected for radial distribution 



P (A" 3 ) 


E/N (K) 


T/N (K) 


P (bar) 


c (m/s) 


0.00352 


-0.040(3) 


4.51(2) 


-0.06(1) 


70(3) 


0.00493 


0.054(3) 


6.84(2) 


0.63(2) 


167(7) 


0.00634 


0.552(6) 


9.68(4) 


2.9(2) 


276(15) 



TABLE III. Results for liquid D^2 at different densities p, 
with backflow correlations included in the model: energy per 
particle (E/N), kinetic energy per particle (T/N), pressure 
(P), and speed of sound (c). Figures in parenthesis are the 
statistical errors. 

functions that we have plotted in Fig |3] for several densi- 
ties, the main peak is shifting to the shorter distances and 
its strength is increasing with the density. The more pro- 
nounced structure weakly emerges for the largest density 
included in the investigation, where a weak indications of 
the second peak formation can be recognized. The inset 
shows g(r) for the density p — 0.00052 A~ 3 and it sug- 
gests that very dilute regime should be represented by a 
monotonic function. 

Similar conclusion about the structure of the system is 
reflected through the results of the static structure func- 
tion S(k) of DJ,i that we report in FigfJJ The reported 
results are the Fourier transforms of the g(r) functions 
for the most of included k. The exception is the region 
of very small k where we have used results obtained di- 
rectly from the calculations. It is evident that the main 
peak increases as the density increases. Also, very dilute 
regime of D4.1 system, shown in the inset, is described 
with monotonic S(k) behaviour. 

B. Bulk liquid DJ, 2 

The bulk D4.2 system is studied in the density range 
from 0.00282 to 0.00634 A~ 3 and in the Table [TTT] we re- 
port for several densities the total and kinetic energy per 
particle. The kinetic energy is calculated as the differ- 
ence between the total energy and the pure estimation of 
the potential energy. In that way the bias coming from 
the choice of the probe wave function used in the simu- 
lations is removed. For all studied densities the obtained 
DMC results are plotted in Fig. [SJ The same analytical 
form ([5]) that we used to fit one part of D^i data, we used 
here to interpolate the equation of state of the liquid DJ,2- 
When the backflow correlations have not been included 
in the model the fitting resulted with the best set of pa- 
rameters e = -0.015(3) K, B = 0.84(4) K, C — 0.52(9) 
K, and po — 0.00367(2) A -3 , and the equation of state 
([5]) is shown as a solid line on top of the DMC data in Fig. 
[5] When the backflow correlations have been included in 
the model, the fitting resulted with the best set of pa- 
rameters e = -0.043(2) K, B = 0.96(2) K, C = 0.56(6) 
K, and po = 0.00381(5) A~ 3 , where the equation of state 
© is shown as a dashed line on top of the DMC data 
in Fig. [51 In both cases the statistical uncertainties for 
the obtained fitting parameters are given as figures in 
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parenthesis. 
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P(A -3 ) 

FIG. 5. Energy per particle of liquid D4.2 without backflow 
correlations (solid triangles) and with backflow correlations 
(solid circles) as a function of the density p. The solid and 
dashed line correspond to the fit to the DMC energies using 
Eq. ([5} • The error bars of the DMC energies are smaller than 
the size of the symbols. 

From the presented results it is evident that with dif- 
fusion Monte Carlo method, in which only the SV cor- 
relations are considered in the probe wave function, it is 
possible to reproduce negative energy per particle for DJ,2 
bulk, i.e., even with simplest model of the correlations 
functions the DMC method reveals DJ,2 as a self-bound 
liquid. We could not obtain negative energy per parti- 
cle using the VMC method and the simplest model of the 
probe wave function in which only SV type of correlations 
are included, similar as authors in Ref.18 commented. In 
the same work Panoff and Clark reported po = 0.004 A -3 
and eo = —0.08(1) K, in case when the probe wave func- 
tion is improved with the optimal Jastrow, triplet and 
backflow correlations. Their reported result for the equi- 
librium energy per particle is lower then our best result, 
e = -0.043(2) K, although they used the VMC method, 
and we used the DMC method. 

The difference between the results can be caused by the 
difference in the improvement of the probe wave function 
model, and because of the including or not including the 
Fermi correction to the VMC results. Namely, from the 
Table U it is evident that number of particles used in sim- 
ulations, and the corresponding Fermi correction for that 
selected number of particles, can dramatically change the 
results. In the density range that was investigated, from 
0.00282 to 0.00634 A" 3 , the value of the Fermi correction 
is similar for 38 and 54 atoms, and is increasing from ~ 
0.05 K to - 0.09 K. For the density p = 0.004 A~ 3 , 
pointed in the RefflU as the equilibrium density, for 38 
and 54 atoms Fermi corrections is ~ 0.07 K. In the Reffl7l 
the same authors cited that they used 54 atoms in VMC 
simulations of the D^ bulk, and it is possible that they 
also used the same number of atoms in the Refll8l for 
D4_2- Since it is not clear if the Fermi corrections has 
been included in the ReflU or not, this should also be 
taken into account as a possible origin of the quantitative 
difference between our results. 

Concerning our results, it is clear that including the 



backflow correlations moves equilibrium density to a 
slightly higher values (from p = 0.00367(2) A~ 3 to 
po = 0.00381(5) A~ 3 ), and expectedly lowers the equi- 
librium energy per particle (from eo = —0.015(3) K to 
e = -0.043(2) K) . The value obtained in our DMC 
calculations for the equilibrium density po = 0.00381(5) 
A -3 expressed in units of er~ 3 is po = 0.191 cr -3 (a = 
3.6892 A) is slightly lower then the one obtained in VMC 
calculations of Panoff and Clark 18 , po = 0.2 cr -3 . 

Since liquid DJ,2 resembles liquid 3 He, it is useful to 
compare the equilibrium densities and energy per par- 
ticle at the equilibrium densities. Our result for the 
equilibrium density of the liquid D4.2, as well as the 
result of the Panoff and Clark, reveals a lower equilib- 
rium density, compared to the one obtained in liquid 
3 He, po = 0.274 a~ 3 (a = 2.556 A). We address this 
difference to the Fermi nature of the systems, because 
smaller difference was reflected in the comparison of the 
equilibrium densities of the Bose liquid TJ. and liquid 
4 He. Namely, the equilibrium density of the liquid T4- in 
terms of a = 3.6892 A is po = 0.375 er~ 3 , and equilib- 
rium density of the liquid 4 He in terms of a = 2.556 A is 
po = 0.365 cr~ 3 . 

In addition, energy per particle at the equilibrium den- 
sity of the liquid 3 He, eo = —2.464(7) K, is significantly 
lower then the eo of the liquid D4.2 (eo ~ —0.1 K). Again, 
we address this huge difference to the very shallow DJ,- 
DJ, interatomic potential, as well as to the Fermi nature 
of the systems. The Fermi nature of the systems reflects 
even more the extreme quantum properties of the sys- 
tems. Namely, in the case of the liquid TJ, we obtained 
eo = —3.656(4) K, which is more comparable to the value 
e = —7.267(13) K obtained for the liquid 4 He. 
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P(A- 3 ) 

FIG. 6. Pressure and speed of sound of D4.2 (solid lines) and 
DJ,3 (dashed lines) as a function of the density. Left (right) 
scale corresponds to pressure (speed of sound). 

The pressure and the speed of sound can be determined 
from its thermodynamic definitions which suppose know- 
ing of the equation of state. Combining the obtained 
equation of state (JSJ) and thermodynamic definitions, we 
derived the pressure of the system as a function of the 
density 

n,> =,*(|), (6) 
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and using obtained pressure we extracted the speed of 
sound as function of the density 

The extracted values for the pressure P and the speed of 
sound c for several investigated densities we included in 
the Table HH In addition, using the Eqs. © and (7J) we 
obtained the functions P(p) and c(p) and we presented 
them in FigJU] The same equations we used to determine 
the spinodal density and the spinodal pressure in liquid 
DJ,2. The speed of sound becomes zero in liquid D4.2 
at density p s = 0.002813 A~ 3 = O.Ulcr" 3 , and at very 
small negative pressure P s = —0.11(1) bar. In terms of 
sigma, the spinodal density in liquid is lower then 
in liquid 3 He ( p s = 0.202er~ 3 ) and even closer to the 
equilibrium density then in liquid 3 He. 

1.2 1 , , , 1 



1.6 




4 8 12 



r(A) 

FIG. 8. Two-body radial distribution functions of the liq- 
uid D^2 for atoms having different spin orientations. From 
bottom to top in the height of the main peak, the results 
correspond to densities 0.00352, 0.00493, and 0.00634 A -3 . 



For the same densities we report in FigO the corre- 
sponding results for S(k). Similar as in D4.1 bulk, the 
main peak increases as the density increases. 



p=0.00352A^ 
p=0.00493 k~i 
p=0.00634 A~ 3 



r(A) 



FIG. 7. Two-body radial distribution functions of the liquid 
DJ,2 for atoms having the same spin orientation. The results 
correspond to densities 0.00352, 0.00493, and 0.00634 A -3 . 

The two-body radial distribution function is also ob- 
tained with pure estimators for liquid DJ.2- For three 
densities in Figfjj we plotted the DMC results for g(r) 
for atoms having the same spin orientation, and in Figj8] 
for atoms having the different spin orientations. In both 
cases the similar behaviour is present, i.e. when p in- 
creases, the structure starts to be more pronounced, that 
can be seen in larger peaks and smaller average first- 
neighbour distances. On the other hand, the compari- 
son between the presented two-body radial distribution 
function reflects the spin-dependent difference which is 
consequence of the Fermi statistics. Since between the 
atoms having the same spin orientation repulsion should 
be more effective, the main peak in FigfJJ at the density 
p = 0.00352 A -3 is practically not possible to localize, 
while for the same density the main peak can be clearly 
localized in Fig|8] Also, due to the effective attraction 
between the atoms having the different spin orientations, 
the main peak at the density p = 0.00634 A" 3 in Figfg] 
is significantly higher (g(r) > 1.2) then the main peak 
at the same density in FigJT] (g{r) < 1.2). The second 
peak emerging can be recognized at p — 0.00634 A~ 3 for 
both cases, the same spin orientation and the different 
spin orientations of the atoms. 





p=0.00352 A ! 




p=0.00493A^ - _ 




p=0.00634A~ J — 











2 

k(A- 1 ) 



FIG. 9. Static structure function of D^2- From bottom to 
top in the height of the main peak, the results correspond to 
densities 0.00352, 0.00493, and 00634 A" 3 . 



C. Bulk liquid D| 3 

The liquid DJ,3 has been also studied in the density 
range from 0.00282 to 0.00634 A" 3 , and in the Table |W] 
we report for several densities the total and kinetic en- 
ergy per particle. For all studied densities the obtained 
DMC results are plotted in FigfTU] and used to interpo- 
late the equation of state given by the analytical form ([5]) . 
When the backflow correlations have not been included 
in the model the fitting resulted with the best set of pa- 
rameters e = -0.181(2) K, B = 0.87(4) K, C = 0.52(7) 
K, and p = 0.00372(2) A" 3 . The equation of state © 
we plotted as a solid line on top of the DMC data in 
FigfTU] When the backflow correlations have been in- 
cluded in the model, the best obtained set of parameters 
is e = -0.229(1) K, B = 1.01(2) K, C = 0.64(5) K, 
and po — 0.00389(2) A -3 , where the equation of state 
is shown as a dashed line on top of the DMC data 
in FigfTU] In both cases the statistical uncertainties for 
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P (A" 3 ) 


E/N (K) 


T/N (K) 


P (bar) 


c (m/s) 


0.00352 


-0.220(2) 


4.34(2) 


-0.08(1) 


65(3) 


0.00493 


-0.143(4) 


6.60(4) 


0.58(2) 


164(5) 


0.00634 


0.327(6) 


9.34(4) 


2.9(1) 


274(11) 



TABLE IV. Results for liquid D4.3 at different densities p, 
with backflow correlations included in the model: energy per 
particle (E/N), kinetic energy per particle (T/N), pressure 
(P), and speed of sound (c). Figures in parenthesis are the 
statistical errors. 

the obtained fitting parameters are given as figures in 
parenthesis. 

The obtained DMC results show that the ground state 
of the DJ,3 is liquid, even when only the SV type of corre- 
lations are used in the probe wave function. PanofF and 
Clark in Refd also reported negative VMC energy per 
particle for DJ,3 , in the case in which only the SV type of 
correlations they used to model the probe wave function, 
indicating in that way the liquid ground state of the sys- 
tem. We reproduced their VMC results using 57 atoms 
in simulations and the same variational parameter b they 
used for the SV correlations, and we noticed that con- 
clusion about the DJ,3 ground state changes if the Fermi 
correction is added to the VMC results. In the density 
range from 0.00282 to 0.00634 A~ 3 , the Fermi correction 
for 57 atoms increases from 0.04 K to 0.07 K. Adding this 
correction to the VMC results does not produce negative 
energy per particle in the above mentioned density range, 
and does not predict the liquid D4.3 ground state when 
only SV type of correlations are included in the model. 

If we compare our DMC results in which the back- 
flow correlations are not included with those in which 
the backflow correlations are included, it is evident that 
including the backflow correlations in the model moves 
equilibrium density to a slightly higher values (from 
po = 0.00372(2) A~ 3 to po = 0.00389(2) A" 3 ), and 
lowers the equilibrium energy (from eo = —0.181(2) K 
to eo = —0.229(1) K), as it was already noticed for DJ.2. 
In addition, the equilibrium density of the DJ,3 is always 
slightly higher then the equilibrium density of the DJ,2- 
The results reported for the DJ,3 equilibrium density in 
RefQJ is slightly higher (p = 0.004 A" 3 ) than one we 
obtained with DMC method. 

Even thought the energy per particle at the equilib- 
rium density eo = —0.229(1) K of the D| 3 is lower then 
the eo = —0.043(2) K of the DJ, 2 , it is still very small neg- 
ative value for the energy per particle at the equilibrium 
density, which defines DJ,3 as very weak self-bound liquid, 
and also very unique in a sense that does not poses its 
helium analogue. The best VMC result for eo = —0.21(1) 
K reported in ReffLsl is similar to our result even though 
we used the DMC method. Including the Fermi correc- 
tion to the VMC results reported in the Reffl8l would 
produce less negative energy per particle since it is evi- 
dent from the Table HI] that the Fermi correction for 57 
atoms is positive and not negligible. 
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FIG. 10. Energy per particle of liquid DJ,3 without backflow 
correlations (solid triangles) and with backflow correlations 
(solid circles) as a function of the density p. The solid and 
dashed line correspond to the fit to the DMC energies using 
Eq. (Jsj . The error bars of the DMC energies are smaller than 
the size of the symbols. 



As in the case of DJ.2, we used the Eqs. (j6]) and © to 
derive the functions P(p) and c(p) for liquid D4.3. The ex- 
tracted values for the pressure P and the speed of sound 
c for several investigated densities are included in the 
Table IIV1 while functions P(p) and c(p) for liquid DJ, 3 
are added in Fig[6] next to the results obtained for liq- 
uid D^2- In liquid D^3 the speed of sound becomes zero 
at density p s = 0.002903 A" 3 = 0.146cr~ 3 , and at very 
small negative pressure P s = —0.12(1) bar. It is evident 
already from the Fig [6] that there is a small difference 
between the spinodal densities in DJ,2 and D4.3 liquids, 
even though the spinodal pressures are rather similar. In 
addition, very similar values of the pressure and speed 
of sound in D4.2 and DJ,3 liquids are revealed from the 
FigH as well as from the Table lED and Table [TV] in the 
remaining region of investigated densities. 




-0.3 1 1 1 1 1 

0.002 0.004 0.006 

P(A- 3 ) 

FIG. 11. Energy per particle of DJ,i (solid squares), DJ,2 
(solid circles) and D4.3 (solid triangles) as a function of the 
density p, with backflow correlations included in the model. 
The presented lines correspond to the fit to the DMC energies 
using Eq. The error bars of the DMC energies are smaller 
than the size of the symbols. 
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IV. CONCLUSIONS 



The ground state properties of the three spin polar- 
ized deuterium species have been accurately determined 
using the DMC method within the fixed-node approxi- 
mation. The accuracy of the DMC method and precise 
knowledge of the DL-DJ. interatomic potential allowed 
for a nearly exact determination of the main proper- 
ties of the systems. The best obtained results for en- 
ergy per particle for all three Dl species are summa- 
rized in FigfTTl The energy ordering for the three Di 
species, close to the equilibrium densities, was found to 
be (E/N) DU > {E/N) DU > (E/N) DU , due to the de- 
generacy, as it was already pointed in previous variational 
descriptions of the systems. The presented equations of 
state indicate that at higher densities the DJ,i should be 
energetically more preferable then D^2- 

Our study confirms previous variational predictions for 
the self bound quantum liquids D4.2 (po = 0.00381(5) 
A" 3 ) and D| 3 (p = 0.00389(2) A" 3 ). The spinodal 
densities are determined also for both liquids. We also 



discuss the ground state of D|i and predict the gas den- 
sity and the pressure at which D^i system could liquefy 
at T=Q, i.e. the conditions at which the system could un- 
dergo a first order gas-liquid phase transition. In order 
to investigate even more interesting properties of the D4.1 
system, we will continue with the diffusion Monte Carlo 
calculations for the solid phase. It will be intriguing to 
find out if D-li is the system that could undergo more 
then one phase transition at T=0. 
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